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Abstract. In this paper, we present a method that enables to solve in parallel 
the Euler-Lagrange system associated with the optimal control of a parabolic 
equation. Our approach is based on an iterative update of a sequence of inter- 
mediate targets and gives rise independent sub-problems that can be solved in 
parallel. Numerical experiments show the efficiency of our method. 

Dans cet article, on presente une methode permettant une parallelisation en 
temps de la resolution des equations d'Euler-Lagrange associees a un probleme 
de controle optimal dans le cas parabolique. Notre approche est basee sur une 
mise a jour iterative de cibles intermediaires et donne lieu a des sous-problemes 
de controle independants. Les resultats numeriques prouvent I'efHcacite de la 
methode. 



1. Introduction: 

In the last decade, time domain decomposition has been exploited to accelerate 
the simulation of systems ruled by time dependent partial differential equations [T] . 
Among others, the parareal algorithm [5] or multi-shooting schemes [1] have shown 
excellent results. In the framework of optimal control, this approach has been used 
to control parabolic systems [2], [I]- In this paper, we introduce a new approach to 
tackle such problems. The strategy we follow is based on the concept of target tra- 
jectory that has been introduced in the case of hyperbolic systems in [6]. Because 
of the irreversibility of parabolic equations, a new definition of this trajectory is 
considered. It enables us to define at each bound of the time sub-domains relevant 
initial conditions and intermediate targets, so that the initial problem is split into 
independent optimization problems. 

We now introduce some notations. Given d € N, we consider the optimal control 
problem associated with a heat equation defined on a compact set fi C K"^ and 
a time interval interval / = [0,r]. The control applies on a subset ftc C ft and 
belongs to the Hilbert space H := L^(/; L^(ric)) whose scalar product and norm 
are denoted by (.,.)•« and ||.||-h respectively. We also denote by ||.|| and ||.||c the 
norms associated with the spaces L^(ri) and L^(ric) respectively. Given a function 
(p defined on / x 17, we denote its restriction to a sub-domain /' x C / x by 

The paper is organized as follows. The optimal control problem is presented in 
Sec. [2]; Sec. [3] is devoted to the description of the time parallelization setting. The 
algorithm is given in Sec. ID Finally, numerical results are presented in Sec. [3 
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2. Optimal control problem 

Given a > 0, T > 0, a target state ytarget G L'^i^), we consider the optimal 
control problem : Find G H 

V* := argmin^g^ J(i;), 

where J is the quadratic cost functional defined by 

(1) Jiv):=lMT)~ytar,etr + ^\\v\\l,. 

The state variable y depends linearly on the control v through the evolution equa- 
tion 

(2) dty - i^Ay = Bv, 

with an initial data j/o G L'^{ft). In this equation A denotes the Laplace operator, 
I' is the diffusion coefficient and B is the injection from L^(i7c) to L'^{fl), so that 
Bv e L^I;L\n)). 

The corresponding optimality system reads 

dty — v/^y — Bv on / x 



(3) 



(4) 



y(0) = yo, 

dtP + J^Ap = on / X 



P{T) = y{T) - ytarget, 

(5) av + B*p^O, 
where B* is the adjoint operator of B. 

Note that for any a > 0, the functional J is strictly convex, so that the system 
(jSHSJ has a unique solution v* . We denote by y* , p* the associated state and adjoint 
state. 

3. Time parallelization setting 

We now aim at solving in parallel the coupled system corresponding to Equations 
(jSHSJ. To do this, we decompose the interval / into subintervals and introduce a 
set of intermediate target states so that the initial problem is replaced by a set of 
smaller independent optimal control problems. The resolution of these problems is 
achieved using an inner loop, whereas the intermediate states are updated by an 
outer loop. 

We start with the definition of the target states. Given a control v € % and its 
corresponding state y and adjoint state p, we define the target trajectory by: 

(6) x(^) = y(^) ~ pI"*^) on / X f2 

In what follows and for the sake of simplicity, we omit the dependence oi y, p and x 
on the control v in the notations. The introduction of this trajectory is motivated 
by the following result. 

Lemma 1. We keep the previous notations. Let r g]0,T[, and the optimal control 
problem: Find w* such that 

w* := argminweHJTiu!), 

where 

Jr{w) := \\\y{T) - x*{r)f + ^\Mh^^^o^r]■L^n.)) 
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with yir) the solution of the Equation We have 

wt = vt, , 

^ I [0.^1 

Proof: Thanks to the uniqueness of the solution of the optimization problem 
associated to Jt , it is sufficient to show that v*^^ ^ is a solution of its optimality 
system. This one is given by Equation (jS]) (with v — Wr) restricted to [0,t] x n 
and 



(7) 



dtp + vAp = on [0, T]xn 
p{t) = y{T) ~ x{t), 



(8) aWr + B*p = 0. 

First, note that y^^ obviously satisfies Equation ([3|) restricted to [0, t]xH with 
V = v*^^ J. It directly follows from the definition of x* (see that: 

p*{T)=y*{r)^X*{r), 
so that p*^^ J satisfies (O . Finally, Equation ([8]) is a consequence of ([U . The result 
follows. □ 
Given > 1, we decompose the interval / — [0, T] into subintervals / ~ U^Jq^/„, 
where /„ = [^„,i„+i], to — < ti < ... < t^-i < t^ ~ T. We also introduce the 
spaces Tin '■= L"^ [In] i^c)) and the corresponding scalar product (.,.)•«„ and 
norm In this framework, given Vn S T-in we define as follows 

(9) < — argmin„^g^_^ J^(u„), 
with 

(10) rM := \\\yn[tn+l) - X{tn+lW + ^Wfu^, 

where x is associated to v through the definition In this functional, the state 
yn is defined by 

dtyn - vAyn = Bvn on /„ x 

These subproblems have the same structure as the original one and are also strictly 
convex. Note also that their definitions depend on the control v through the target 
trajectory, hence the notation J^. 

The optimality system associated with these minimization problems are given by 
Equation pT|) and 



(11) 



(12) 



dtPn + vApn = on /„ X 

Pn{tn+l) = y{tn+l) - X{tn+l), 



(13) avn + B*pn = 0. 

Lemma 2. We keep the previous notations. Denote by x* the target trajectory 
defined by Equation ^ with y = y* and p = p* and by y^,p^,v^ the solutions of 
Equations Ul\]13\) associated with v* . One has: 
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The proof of this result foUows the hnes of Lemma [Hand is left as an exercise to 
reader. 



4. Algorithm 

We are now in the position to propose a time parallelized procedure to solve 
Equations ([3])-([5]). Consider an initial control and assume that, at step k, a 
control v'' is known. The computation of is achieved as follows: 

(1) Compute y'', -p^ and the associated target trajectory ^ according to Equa- 
tions (121), (01) and © respectively. 

(2) For n = 0, ...,iV — 1, solve the sub-problems ^ in parallel and denote by 

the corresponding solutions. 

(3) Define v^^^ as the concatenation of the sequence ('5^'*'^)ri=o,...,Af-i- 

(4) Update the control variable by 

(14) w'^+i =t;'=-f0fe(^'=+i 

where the value Qk is chosen to minimizes J{v^ -I- dk{v''~^^ ~ v'^))- 
We have not detailed Step [2] as we rather aim at presenting a structure for a general 
approach. However, because of the strict convexity of the problems we consider, 
a small number of conjugate gradient method steps can be used to achieve the 
resolution of these steps. 

5. Numerical experiments 

In this section, we test the efficiency of our method. We consider a 2D example, 
where = [0,1] x [0,1] and ilc = [j: |] x [j; |] ■ The parameters related to 
our control problem are T = 6.4, a = 10"^ and u = 10^^. The time interval is 
discretized using a uniform step St — 10"'^, and an Implict-Euler solver is used 
to approximate the solution of Equations (jSHl]). For the space discretization, we 
use Pi finite elements. Our implementation make use of the freeware FreeFem [8] 
and the parallelization is achieved thanks to the Message Passing Interface library. 
The independent optimization procedures required in Step [2] are simply carried out 
using one step of an optimal step gradient method. The results are presented in 
Figure [1] 

In the first plot, we consider the evolution of the cost functional values with 
respect to the iterations and do not take into account the parallelization. The 
result reveals that our algorithm significantly accelerates the optimization process. 
This outcome may indicates that the splitting introduced in our approach acts as 
a preconditionner during the numerical optimization. This will be the purpose of 
some further investigation [9], in the same spirit as in [?]. 

In a second plot, we represent the evolution of the cost functional values with 
respect to the number of matrix- vector product. Parallel computations that are 
done in Step [2] are only counted once. When comparing with a standard optimal 
gradient step method, we observe speed-up approximativcly equal to 3. 
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Figure 1. Functional values evolution, with respect to the num- 
ber of iterations (top) and multiplications (bottom). 
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